Cast Shadows and Linear Subspaces for Object 

Recognition 

Cross Reference to Related Applications 

[0001] This Application is a Continuation -in -Part of U.S. 
Patent Application Serial Number 09/967,206 (Attorney 
Docket No. NECI1110, 14919), filed 28 September 2001, 
entitled "Broadened- Specular Reflection and Linear 
Subspaces for Object Recognition", the disclosure of which 
is hereby incorporated by reference for all purposes. 

Background of the Invention 
Field of Invention 

[0002] The present invention relates generally to machine 
vision, and more particularly, to a method of modeling the 
shadows cast by various features of an object onto its own 
surface under certain lighting conditions, to allow a 
machine to recognize the image of that object when 
illuminated by a wide variety of lighting conditions. 

Description of Related Art 

[0003] It is clear that changes in lighting can have a 
large effect on how an object, such as a person, looks. 
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Understanding this is essential to building effective 
object recognition systems. 

[0004] U.S. Patent Application Serial Number 09/705,507 
(Attorney Docket No. NECI 1093, 13919), filed 3 November 
2000, entitled xx Lambertian Reflectance and Linear 
Subspaces" , the disclosure of which is hereby incorporated 
by reference, considered the relationship between the 
function that describes the lighting intensity, and the 
reflectance function that describes how much light an 
object reflects as a function of its surface normal, under 
a given lighting condition. Representing these functions 
as spherical harmonics, it was shown that for Lambertian 
reflectance the mapping from lighting to reflectance is a 
convolution with a nearly perfect low-pass filter. High- 
frequency components of the lighting hardly affect the 
reflectance function. Therefore, nearly all Lambertian 
reflectance functions could be modeled as some linear 
combination of nine spherical harmonic components. 

[0005] These low-dimensional approximations have two key 
advantages. First, they are described with few parameters, 
thus optimization is simpler. Moreover, the images 
produced under low- frequency lighting are inherently immune 
to small changes in orientation. A small rotation in the 
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surface normal is equivalent to a small rotation in the 
light. However, low-frequency light does not appreciably 
change with a small rotation. The analysis shows this, in 
that an image produced by an n th order harmonic light 
component is computed by taking an n th degree polynomial of 
the components of the surface normal. 

[0006] U.S. Patent Application Serial Number 09/967,206, 
incorporated by reference above, expanded the spherical 
harmonic approach to consider both the effects of broadened 
specular reflection as well as Lambertian reflection. 
Thus, the reflectance function that describes how much 
light an object reflects as a function of its surface 
normal, under a given lighting condition was made more 
realistic . 

[0007] That work has considered the effects of attached 
shadows, which arise locally, when some portion of the 
object is facing away from the light. However, it has not 
addressed the problems that arise when one part of an 
object casts a shadow on a different part of the object. 
This is important, because many objects of interest are 
characterized by features whose prominence produces cast 
shadows over a significant portion of the surrounding 
surface, especially at low angles of illumination (high 
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angles of incidence relative to the surface normal) . 
Indeed, no natural object is convex, and relatively few 
manufactured objects are, and even fewer when seen in 
combination. Thus, cast shadows can be expected to play an 
important role in reflection, and, therefore must be 
included in harmonic treatments. 



Brief Summary of the Invention 

[0008] The goal of the present invention is to allow a 
machine, such as a computer, to match an input image of an 
object, such as a human face, to an object model stored in 
a database associated with the machine, and to perform this 
matching with a high degree of accuracy. 

[0009] The steps of the recognition method of the present 
invention include providing a database of three-dimensional 
object models, providing an input image, positioning each 
three dimensional model relative to the input image, 
determining an optimal rendered image of the model that is 
most similar to the input image, computing a measure of 
similarity between the input image and the optimal rendered 
image of each model, and selecting the object whose optimal 
rendered image is most similar to the input image. The 
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step of determining an optimal rendered image includes 
deriving a reflectance function for each object that 
includes the effects of cast shadows, and optimizing that 
reflectance function to most closely match the input image. 

[0010] The method of the present invention includes 
rendering the images produced by an object model when 
illuminated only by each of a plurality of spherical 
harmonic components of light. The summation of these 
components defines the reflectance function for the object 
model. The derived reflectance function takes into account 
the intensity of light components incident on the object 
model, as well as the portion of incident light which does 
not reach a given point because that light originates below 
some local horizon, that point therefore being in a region 
of cast shadow. The horizon is defined in terms of the 
polar angle at a given point, and can be a function of the 
azimuthal direction of the incident light. 

[0011] Among the advantages of this approach is the ease of 
application and computational efficiency over serial 
methods of dealing with cast shadows, such serial methods 
being of the type including ray tracing, for example. 
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Detailed Description of the Invention 
Deriving the Reflectance Function 

[0012] To derive a reflection function that accounts for 
some obstruction between the surface point of interest and 
the light source, consider a small surface region of an 
object model with surface normal n. Let \i = cos9, where 9 
is the polar angle of a ray of light striking the region 
relative to n , Let 0 be the azimuth of the ray, also 

relative to n. Now if the object is convex, all rays with 
O<0 <k/1(/u>0) and 0<(f><2x will illuminate the surface. 
However, in general, only rays whose polar angles satisfy 
ju(0)<ju<1, O<0<0(</>) will contribute to the illumination, 
and hence to the reflection. Here 0{(j>) is just the angle 
(relative to the normal) that outlines the local horizon as 
one sweeps through the azimuth. In the usual treatments 
one takes e($) = 7zll , for if 0(§)ktiI1 , so little light 
intensity per unit area is available that, unless the 
incoming light is low and highly directional, little error 
is made. However, for 0{<j>) much less than nil, some 
compensation must be made for the resulting shadows. 
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[0013] Three models of cast shadows are considered in the 
context of the harmonic -decomposition approach. The most 
general model takes 0{<j)) and determines the surface 

intensity for each L lm lighting component, where 

L(n i ) = J Z lm L lm Y kn (n i ) is the illumination on the object. This 

involves considerable calculation, but will be necessary in 
many realistic examples. 

[0014] The simplest model, the crater model, is to assume 
that ju(<fi) = ju 0 . In other words, that the region of interest 

sits in a depression with symmetry about the azimuth. 
While it is not clear, without more computation, what 
effective value to choose for ju 0 given ju(<f>) , for ju 0 >0 

clearly the reduction in incident intensity and shift to 
higher harmonics can at least be simulated. These features 
are characteristic of shadows and are impossible to 
incorporate treating objects as convex. Sketching the 
derivation of the reflection in this model following the 
lines and notation of Application Serial Number 09/967,206, 
incorporated by reference above, one can then consider a 
somewhat more general model, and finally proceed to the 
most general case. 
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[0015] If r cscx (h) is the reflection in the presence of cast- 
shadows in the crater approximation for a surface normal of 
n , then 



where L(n s ) is the incident light intensity in the -n. 
direction and ju 0 defines the effective local horizon of the 
crater. In the 09/967,206 Application ju 0 was taken as 

zero. Also, u denotes the step function that is one for 
positive values and zero otherwise. Therefore, the step 
function including ju as shown in the above equation 
analytically excludes light originating from any point 
below the local horizon at the point of interest from 
contributing to the reflectance. 

[0016] Focusing on the Lambertian factor one can find 



where Y lm (n) is the lm spherical harmonic function of the 
direction n, and where 
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[0017] Since 



Si 



^(^^) = ((2/+l)/4^) 1/2 i>(//) 
^=((2/+l)^) 1/2 f 1 rf/vilJC^) 



evaluating one can find 



J/2/1 .,3w>l/2 



u ^=^ 2 (l-^)/3' 

5 * -rr2/+n.Ty /2 f (/ " 2)!( " 1)(/ " 2)/2 is' rMWJJ^^-^ 

1 ^-tt^+W < 2 '({-l)!(£^)! W 2'Jt!</-Jt)K2it-C/-2)>! ^° } 



for />2, / even, and 



/tf mw)/2 2 l k\(2k-(I-2))\ ^ 



for />3, / odd. For // 0 =0 the previous results are 
recovered. For ju 0 >0 , there is no first-order (in ju 0 ) 
change in k^ . Second-order changes occur in k M0 and k^ for 
2 even: k M2 increases as -((5;r) 1/2 14)$ for ju 0 small. 
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[0018] In the foregoing, the k Ml coefficients are shown for 

the harmonic expansion of the light reflected from a 
Lambertian surface in the presence of cast shadows in the 
crater approximation, where the elevation of the local 
horizon ju(<f>) = cos(6>($)) is taken to be independent of azimuth, 

= m a more general model, the full azimuthal 

dependence ju(<f>) is included, but the Lambertian factor is 
assumed to have the same form as in the crater model. 
While this restricts the actual dependence of the 
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reflection on the surface normal, it provides an indication 
of the role played by higher harmonics. 

[0019] Replacing ju 0 with ju(<f>) , the solution is 
straightforward to obtain from the crater result. Consider 
k^i and treat it as a functional of ju(</>) rather than as a 

function of ju 0 . Thus 

where ju = cosO . For fi{(j>) = ju 0 , is equal to 2tc times a 
polynomial of degree (/ + 2) in ju 0 . For general ju((f>) , replace 
ful in the expression for by 

Vl-*\ 2 *d<l>lf{<l>)l2n 

[0020] If l L is the value of the index / of the largest 
harmonic of interest, we need only consider n such that 
2<n<(l L +2):l L =1 for the four -dimensional model (i.e. the 0 th 
and 1 st order spherical harmonic components); l L =2 for the 
nine -dimensional model (i.e., the 0 th , 1 st and 2 nd order 
spherical harmonic components) . Note that the n=0 case is 
trivial, and that the n=l case does not enter. 



-11- 



G:\Nec\l 196U49 1 9z\spec\149 19Z.spec.doc 



[0021] Even though the table given for k pl assumes ju(</>) = ju 0 , 

it is actually useful in the more general case. To see 
this let 



[0022] Then, since ju tt > jf for n>\ , note that the actual 
values of k Ml for ju($) are indicative of the values given for 

k Ml for //(^) = // 0 , where ju<ju 0 . While admittedly rough, it is 
at once apparent that the dependence of ju on </> enhances 

the harmonic content of the image. Also, since ju> ju 2 > //... 
for most cases of interest, only relatively low values of n 
(powers of ju ) will play a significant role in the £ , . 

[0023] To include the full effect of cast shadows in the 
harmonic representation, one must fully express the 
expansion of the Lambertian factor. In both the crater 
model and the more general model an assumption was made to 
ignore the azimuthal dependence in the unit step function 
//(•) , the (j> dependence of ju(<f>) . While satisfactory in those 
cases where ju($) is nearly constant, in non-convex regions 
its variation with azimuth can be significant. 
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[0024] Though tricky to evaluate, the solution can be 
written down by analogy with the more general of the two 
cases evaluated above* If r cs (n) is the reflection in the 

presence of cast shadows due to a local horizon of 
ju(<f>) = cos(0(<f>)) that is due to an elevation dependent on 
azimuth <j) , then we can write much as above that 



[0025] Here D l m , m (afiy) is the rotation matrix associated with 
spherical harmonics, and (fi,a) are the polar coordinates of 
n(y-0). Angles {apy) are the Euler angles of the effective 
rotation that takes the z -axis from the viewing direction 
to the direction of the surface normal n . 

[0026] The evaluation of k lm is complicated by the azimuthal 
dependence of the elevation ju(<f>) . One simplification that 
retains much of the azimuthal or phase coupling introduced 
by cast shadows is to assume that we can approximate ju(<f>) by 

two values, ju x for fa<$<$2 and ju 2 for </> 2 <$<$ } +2x . Then 
the integrals over d$ and dju separate, and one obtains 
factors of sin{m{(j> 2 -<l> x )l2)lm entering the strength of the 
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contribution to from the fi($) = fa and ju($) = <f> 2 regions. 
Phases pn l(<f> 2 -<i>x) , p odd, yield the largest k lm . Thus for 
m = 0 

and for m>0 or m<0 

k lm * J*sin{n4ix /2)/(m/2)-r 2 djujuY lm (0) 

where ? = (^+^)/2,^ =tt«*)A(«) = c* l> 7 fa (^) f which is 
independent of ^ . 

[0027] Another simplification arises when the local horizon 
is axially symmetric about an axis m other than the surface 
normal n; we can call this the megaphone approximation and 
note its potential utility when only n and m , but not the 
megaphone aperture ju m , vary with surface position. We 
write 

r CSM (n;m^ m ) = Jrfn^Xn, n v0)w(n, -m»// m ) 

(the expression Aw B is equal to the greater of (A, B) ) and 
express the three factors in the integrand as 
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L(n f ) = S /w AA(A f ) 
(n,..nvO) = S /in)i ^^(n)^(n ; .) 

where k, = Cd{i{iY l0 (ju) and # = f 1 dfiY l0 ^) . By this point the 

analogy with diffuse- specular reflection should be clear. 
The term p l has replaced kf , but otherwise the solution to 
that problem can be taken over intact. Thus one finds 

r CSM (n; m, //J = £ /w A^« 2 m 2 J (W ; ™i , ™ 2 > m)k hPl Y lm (m)7^ (n) 

where 

/(/ 1 / 2 /;m 1 ,m 25 m) = Jd^ft)!^ (4,)!^ ft) 

= ((21, +l)(2/ 2 =l)/4^(2/ + l)) 1/2 C(/ 1 / 2 /;m 1 m 2 m)C(/ 1 / 2 /;0 

where C{l x lJ\m x m 2 m) are the Clebsch-Gordon (C-G) coefficients. 
When applicable, this approach is significantly more 
computable than the general solution above, increasing the 
speed of computation and/or maximizing the use of 
processing resources. 

[0028] The foregoing expressions for the modification of 
the reflected intensity in the presence of cast shadows 
requires that the horizon fi{<j>) be determined for each 
surface position of interest. This requires two steps. 
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First, the surface must be rotated so that the normal n at 
this position is parallel to the z axis. Then the maximum 
elevation ju($) must be found in each direction of interest 
<j> . ( ju($) = cos9{(j>) where 0(<t>)is the minimum polar angle from 
the vertical in the direction <j> . ) The assumption is made 
that a set of surface positions {r} = { (r x ,r y/ r z ) } is given in 
some known frame of reference. 

[0029] The first task is to rotate each surface position r 
into r' such that the normal n at the position of interest 
r 0 becomes parallel to z at r 0 ' . This is accomplished by 
letting r'=Rr, where R is given by 

R = zn r + (fl x Z )(fl x z) r + (z x (fi x z))(n x (fl x z)) T 

[0030] Note especially that Rn = z. The three terms in R 
operate on the component of r parallel to n , to nxz 7 and 
to nx^xz), respectively. With this preliminary step out 
of the way, one can determine the horizon //(^)for position 

[0031] It is now simplest to pass to a set of variables 
relative to r<;=Rr 0 : 
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r ' - r o = (K - r L > r 'y - r o y > K ~ r L ) = (d x ,d y , h) 
d x = (d 2 +d 2 ) l,2 cos0 
d y ={d 2 x +d 2 y ) m sin<l> 

[0032] The (d x ,d y ) give the position of the height h(d x ,d y ) of 

the surface above the plane through r 0 normal to n. It 
follows that ju(<f>) is given by one of the following: 



max(l + d 2 x /(hcos<f>) 2 y h2 
d x 

mca(l + d 2 y /(hsin^) 2 y V2 
0(zero) 

where only d x or d y where h(d x ,d x tan<f>) or h(d y cot<f>,d y ) , 

respectively exceed zero enter into the maximum. For h<0 
for all positions along ^ 5 // = 0. Clearly, one will choose to 
use the first expression for -7tlA<(t><7ilA and for 
2>7ilA<(j><57tlA , and the second for ;r/4<^<3;r/4 and for 
5^/4<^<7^/4. 

[0033] As an example, see the following illustration 
computing ju((f>) for a ridge of isosceles-triangular cross 
section of height h > 0 and base b, and unspecified length 
- somewhat like the nose on a face. At a distance d>b/2 to 
the left of the centerline of the base of the ridge, taken 
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in the direction # = ±n/2 9 ft(f) = (l + (d Ihcos^f)' 112 for -n I2<<j><x 12 
and //(^) = 0 for n I2<<f> <?>n 12 . On the left side of the ridge 
itself ju(<t>) = (l + (b/2h) 2 y m for ;r/2 <<£ < 3;r/2 and //(^) = 0 for 
-7il2<<t><7il2 . The latter is an example of /u(<f>) constant 
over ranges of <j) , whose solution was given above* The 
former can be approximated over appropriate ranges of <j> , or 
integrated over <j> if greater precision is required. Very 
similar solutions apply for a trapezoidal ridge, with 
ju(</>) = 0 along the top of the ridge. 

[0034] The crater model is sufficiently simple to introduce 
most of the effects due to cast shadows. Odd orders 
(1=3,5,7, ...), which play no role in the case of convex, 
Lambertian surfaces, become appreciable as the crater 
deepens. The strength of orders />4 oscilate in value as 
the fit for the various modes waxes and wanes with 
increasing ju 0 . From the general solution, note the 

couplings between the azimuthal frequencies in the horizons 
and those in the resulting expansions. The actual 
calculation of the horizon was indicated only briefly, as 
it has probably been more fully developed in the graphics 
community. 
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[0035] The rotation matrices used in the method are 
discussed by Rose, M.E., Elementary Theory of Angular 
Momentum. They enable the spherical harmonics in a rotated 
coordinate system to be expressed in terms of those in the 
reference system. Thus, 



f\ f\l/2 



D , (aB ) _v (-iy((/ + m)!(/-m)!(/ + mQ!(/-m')l) 
* PT ' t\{l-t-m')\(! + m-t)\(m' + t-m)\ 



■e im ' a e™ 7 (cos(0 / 2)) 



2l+m-m'-2t 



{-sin{f3l2)) 



m'-m+2t 



[0036] In the summation over t, t is such that the four 
factorials in the denominator are well-defined. Thus each 
of t>Q,l-m'>t,l + m>t,t>m--iri , must be satisfied for each 
possible t. Also, l>\m\ and l>\m'\. 

[0037] Clearly D* 0 (a/3y) = 1 . 

[0038] Tables for 1=1 are commonplace: 
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where C = (1 + cos/3)/ 2, S = (l- sin/3) / 2, c = cos/3, s = sin/3 . Values for 
1 = 2 are also needed, which have been worked out in the 
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following table. The e** a and e imy factors are omitted the 
clarity of the table although, of course, these must be 
included in any solution. 
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[0039] Further, a number of indefinite integrals enter in 
connection with the integration over the elevation (ju = cos0) 

of the spherical harmonics. The following table should be 
sel f -explanatory : 
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Application of the Reflectance Function 

[0040] The method of deriving the reflectance function 
described above is used in the object recognition process. 
To match an input image with an object model in a database, 
the object model must first be positioned to match the 
position of the input image. There are standard methods of 
performing this function known in the art, though a 
preferred method is disclosed in U.S. Patent Application 
Serial Number 09/538,204 (Attorney Docket No. NECI 1083, 
13414), filed 30 March 2000, entitled "Method for Matching 
a Two Dimensional Image to One of a Plurality of Three 
Dimensional Candidate Models in a Database 77 , the disclosure 
of which is hereby incorporated by reference. 
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[0041] Once positioned, images of the model are rendered 
under each of one or more spherical harmonic components of 
light, based upon the reflectance function described above. 
These component rendered images can be linearly combined to 
arrive at an optimal rendered image of the model. An 
optimal rendered image means a rendered image that most 
closely approximates the input image. The comparison of 
the various rendered images to the input image is made by 
calculating the Euclidian Distance between the images, in a 
manner that is known in the art. The optimization is 
performed by a least squares fit technique between the 
vector that is the reflectance function and the matrix that 
represents the input image, also known in the art. 

[0042] Other steps follow the optimization in the 
recognition process. One is computing some measure of 
similarity between the input image and the optimal rendered 
image for each model. This is accomplished by measuring 
the Euclidian Distance between each optimal image and the 
input image. The object model whose optimal rendered image 
is nearest the input image, assuming some threshold value 
of similarity must be reached, is a match to the input 
image . 



-22- 



G:\Nec\l 196\14919z\spec\l4919Z.spec.doc 



[0043] While the reflectance function of the present 
invention is described in application to the field of 
machine vision, it is not limited to that field. There are 
numerous parallels between this work and advances in the 
field of computer graphics, and our invention has 
applicability there as well. The field of computer 
graphics is continually striving to produce more realistic 
rendered images based on object models that exist only in 
the computer. The reflectance function of the present 
invention advances that work as well. 

[0044] The invention has been described herein with 
reference to a particular exemplary embodiment. Certain 
alterations and modifications may be apparent to those 
skilled in the art, without departing from the scope of the 
invention. The exemplary embodiment is not meant to be 
limiting on the scope of the invention, which is defined by 
the appended claims. 
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